Image Warping with Scattered

Image warping has many applications in art as well as in image processing. Usually, displacements are computed with mathematical functions or by transformations of a triangulation of control points. Here, di erent approaches based on scattered data interpolation methods are presented. These methods provide smooth deformations with easily controllable behavior. The usefulness and performance of some selected classes of scattered data interpolation methods in this context is analyzed. Zusammenfassung Bildverzerrungen werden sowohl in k unstlerischen Anwendungen als auch in der wissenschaftlichen Bildverarbeitung ben otigt.  Ublicherweise wird die Verschiebung der einzelnen Bildpunkte mit analytischen mathematischen Abbildungsfunktionen oder aus einer Triangulation der Kontrollpunkte gewonnen. In der vorliegenden Arbeit wird ein anderer Ansatz prasentiert, der auf Scattered{Data{Interpolationsmethoden basiert. Dieser Ansatz erzeugt glatte Deformationen mit einfach kontrollierbarem Verhalten. Die Brauchbarkeit und E zienz verschiedener Scattered{Data{Methoden in diesem Zusammenhang wird verglichen. Image Warping with Scattered Data Interpolation Methods Detlef Ruprecht and Heinrich M uller 1 Digital Image Warping The term \image warping" describes methods to deform images to arbitrary shapes. Digital image warping has received a lot of interest in recent years. Much of this interest is due to its large range of applications. Some of them are in \artistic" areas, like computer animation. Here, image warping is a basis for two{dimensional morphing, the smooth transition between keyframe images. Other \artistic" applications include facial animation and free{form deformation of images. Another important area closely related to image warping is texture mapping, used for photorealistic rendering of surfaces in computer graphics [12]. Other applications are in scienti c image processing. In image data acquisition, the acquisition method often causes deformations of the acquired image, e.g. through lens distortion. In satellite imaging, distortions are caused by surface curvature and oblique viewing angles. In ultrasonic imaging, distortions are caused by the varying speed of sound in di erent materials. These deformations have to be recti ed to obtain correct coordinates. This process is called \registering". In registration applications, correspondence points between the deformed image and a reference image might be chosen manually or automatically. Considerable work has been done on automatically establishing correspondences, see e.g. [1, 9, 13, 20]. Another application related to registration is the normalization of samples to a standardized form to allow comparisons between the samples to be independent of size and slight form variations. A medical application is the construction of voxel models from planar slices. The distance between slices is usually much larger than the spatial 1 resolution within each slice. For rendering and for surface reconstruction, e.g. with the marching cube algorithm, however, it is preferable to have voxels with edges of equal length. To achieve this, some interpolation between slices is necessary. This is commonly done with linear interpolation between slices. This method is fast, but the result is not completely satisfactory. By deforming one slice towards the next, a better interpolation can be achieved. The process of warping an image can be divided in two steps: The desired displacements of all pixels in the source image have to be computed. The image has to be resampled to create the output image. The rst step is usually done either by applying a global analytic mapping function to the image pixel positions, or by using a set of control points which specify the displacements of some points in the initial image. From these control points, a triangulation can be obtained and new locations computed by transforming each of the triangles [7, 8]. For the second step, a number of algorithms have been developed in recent years. A simple{minded approach would be to take the coordinates of four neighboring pixels, compute their deformed coordinates with the mapping function from above, and perform a bilinear ll for all pixels within the obtained rectangle ([5], p. 249). Much e ort has been spent on more e cient algorithms for the resampling step. In particular, by separating the resampling step into horizontal and vertical displacements, e ciency can be improved considerably [19, 22]. For a detailed introduction to image warping methods, see [21]. This article focuses on the rst step of image warping, the construction of a suitable mapping function. We will call this the \deformation" stage. It will be pointed out that triangulation based deformation is just a special case of scattered data interpolation methods applied to deformation. Triangulation based methods have the disadvantage of having non{continuous derivatives and not being de ned outside the convex hull of the control points. For aesthetically pleasing results, a large number of triangles is needed. Other deformation methods are described which overcome these limitations. Their advantages and disadvantages are compared. For the resampling stage, routines from the Khoros visualization system [16] are employed. These perform a bilinear ll as outlined above, 2 and are therefore not very e cient. However, they are readily available in a user{friendly interactive image processing environment. As an example, g. 1b shows some subtle and some not so subtle deformations of the classic \Lenna" image, with the control point assignments shown in g. 1a. a) b) Figure 1: Deformations of a woman. a) A well known image from image processing papers, with control point assignments. b) The image deformed with multiquadrics. The remainder of this paper is organized as follows. First, the relation between the problem of image deformation and the problem of scattered data interpolation methods is established. Several scattered data interpolation methods along with their application to image warping and discussion of their suitability follow. After that, the performance of the various methods is compared. 2 Image Deformation and Scattered Data Interpolation The problem of image deformation can be formulated as follows: Image Deformation Problem: Input: n pairs (pi;qi) of control points, pi;qi 2 IR2, i = 1; . . . ; n. Output: An at least continuous function f : IR2 ! IR2 with f (pi) = qi, i = 1; . . . ; n. 3 Besides continuity, other properties of the function f are desirable which concern the visual pleasantness and which are hard to describe in an exact manner. These properties will become recognizable in the following discussion of concrete propositions. The problem of scattered data interpolation is to nd a real valued multivariate function interpolating a nite set of irregularly located data points. For bivariate functions, this can be formulated as: Scattered Data Interpolation Problem: Input: n data points (xi; yi), xi 2 IR2, yi 2 IR, i = 1; . . . ; n. Output: An at least continuous function f : IR2 ! IR interpolating the given data points, i.e. f(xi) = yi, i = 1; . . . ; n. As has been pointed out before [7, 19], the mapping functions can be split into their components fj : IR2 ! IR, j = 1; 2, and each component treated separately. These components can then be interpolated with scattered data interpolation methods. The data points to be used are (pi; qi;j), qi;j the j-th component of qi. 3 Triangulation Based Methods Scattered data interpolation through a triangulation of the data points is a classic approach in scienti c visualization. This method consists of rst dissecting the de nition space into a suitable set of triangles with the given data points being the corners of the triangles. Then, each of the triangles is interpolated independently. Several criteria for an \optimal" triangulation are known [15]. One optimal triangulation, which is based on a Dirichlet tessellation of the data set, is called \Delaunay triangulation". It can be computed with a divide{ and{conquer algorithm of algorithmical complexity O(n logn) where n is the number of data points. If the pixels lie on a regular grid, as is usually the case unless a previous deformation has already caused irregular pixel locations, the correspondence between triangles and pixels is trivially known in constant time, and thus each pixel can be mapped to its new location in constant time as well. As computing the new pixel location is the governing process in deformation, this gives an algorithmical complexity of O(N ) for triangulation based deformation, where N is the number of points to be mapped. 4 a) b) c) d) Figure 2: Deformation with triangulation based methods. a) Delaunay triangulation of the initial positions of all control points. b) Deformation of the triangulation. c) Deformed test grid with linear patches within the triangles. d) Deformed image with linear patches. A continuous, although not smooth, interpolation can be obtained by linear interpolation using barycentric coordinates within each triangle. This method has been applied to image warping in [7]. The visual appearance of the result is quite acceptable if the deformations are small and if enough data points are provided so that changes of the transform coe cients between neighboring triangles are small. A smooth deformation can be obtained by using non{linear patches within the triangles. 5 In [8], a method is described which uses cubic functions obtained from a Clough{Tocher subtriangulation. A problem common to all triangulation based methods for image warping is that foldover can easily occur. The term \foldover" describes the occurence of overlapping deformations, i.e. several non{adjacent pixels in the input image are mapped to the same pixel in the output image. With triangulation based methods, this happens if the orientation of the corner points changes for any of the triangles, i.e. the triangle is ipped over by the deformation. Fig. 2 shows an example of deformation by triangulation with linear interpolation within the triangles. The Delaunay triangulation is shown along with the images. Note the edges in the deformed image where control points next to each other have been deformed dissimilarly. In the lower left corner, an example of foldover is observable. 4 Inverse Distance Weighted Interpolation Methods Inverse distance weighted interpolation methods have been originally proposed by D. Shepard [18], and improved by a number of other authors, among them [3, 4, 6, 14]. For each data value i, a local interpolant fi(x) : IR2 ! IR with fi(xi) = yi; i = 1; . . . ; n is determined. The interpolating function is a weighted average of these local interpolants, with weights dependent on the distance of the observed point from the given data points,f(x) = n X i=1wi(x)fi(x); fi(xi) = yi; i = 1; . . . ; n: (1) wi : IR2 ! IR is the weight function and yi 2 IR the data value at the data point xi 2 IR2. The weight functions must satisfy the conditions wi(x) 0; i = 1; . . . ; n; n X i=1wi(x) = 1; and (2) wi(xj) = ij, i.e.wi(xi) = 1 andwi(xj) = 0; j 6= i; i; j = 1; :::; n: (3) These conditions guarantee the property of interpolation. 6 In [18], the following weight function was proposed: wi(x) = i(x) Pnj=1 j(x) ;with (4) i(x) = 1 (di(x)) ; di(x) the distance between x and xi: (5) The smoothness is determined by the exponent . > 1 assures continuity of the derivatives, with the value of the rst derivative vanishing at the data points. The original approach in [18] used just the value yi at the data points for the local interpolant, i.e. fi(x) = yi. This causes all points away from the data points to be moved towards the average of all yi. While this might be tolerable for interpolation, it is unusable for image deformation, as it causes all points to be warped towards the center of the image, as demonstrated in g. 3a. For this reason, linear or quadratic polynomials are normally used as local interpolants. To determine the local interpolants, estimations of derivative values at the data points are needed. In [4], it has been proposed to compute the local interpolant fi(x) by minimizing the squared error of the mapping of other nearby control points xj with fi. In our experiments, we utilize all other control points, weighted with a coe cient wij to attenuate the in uence of distant points. The corresponding error function Ei(f) is Ei(f) = n X j=1;j 6=iwij kf(xj) yjk2: (6) The optimum local interpolant is obtained by minimizing Ei(f). The coe cients wij should depend on the distance between xi and xj in a manner similar to wi(x), but it is not possible to use wi because of (3). A simple possibility is to use wij = i(xj). This approach has been used in our experiments with local interpolants. The application of inverse distance weighted interpolation to image warping gives f (p) = n X i=1wi(p)fi(p); fi(pi) = qi; i = 1; . . . ; n: (7) fi : IR2 ! IR2 are the local interpolants. They can be rewritten as fi(p) = qi +Di(p pi), with Di : IR2 ! IR2; Di(0) = 0: (8) 7 A simple possibility is to choose theDi as linear transformations. Then Di can be represented by a 2 2{matrix (di;kl), k; l = 1; 2. The error function corresponding to a linear transformation D is Ei(D) = n X j=1;j 6=iwij kqi + 0@ d11 d12 d21 d22 1A (pj pi) qjk2 (9) = N X j=1;j 6=iwij[(d11(pj;1 pi;1) + d12(pj;2 pi;2) + qi;1 qj;1)2 + +(d21(pj;1 pi;1) + d22(pj;1 pi;1) + qi;2 qj;2)2]: (10) The minimum of the error function is obtained with the partial derivatives with respect to the dkl; k; l = 1; 2. These are linear functions of the dkl. Thus, we obtain a system of four linear equations with four unknowns, which can easily be solved for the dkl. For a detailed description of deformation with inverse distance weighted methods, see [17]. With a linear local interpolant, the deformation has linear precision, i.e. if the requested displacements of the control points can be satis ed by a linear transformation, this is what will be obtained. An advantage of inverse distance weighted methods is the large number of ways they can be tailored to speci c needs. The range of in uence of the control points can be limited by weight functions which vanish outside a distance R [4]. The interpolation condition can be relaxed by adding a damping term to the weight function, or strict adherence to the local interpolant can be enforced within a distance R around control points [17]. Many parameters can be tweaked to control the behavior of the deformation. Fig. 3 shows examples of image deformations through inverse distance weighted methods with linear local interpolants and various values for . The results are clearly not too satisfactory. The main reason is that the transition between the range of in uence of the control points is rather uneven. This causes unwanted bends in the deformation, as can be seen in the shoulder and arm area of the picture. Experiments with other weight functions and other exponents have not led to signi cant improvements. The best results were obtained with locally bounded weight functions [4], g. 3d, but even this image shows some unevenness, and the quality depends strongly on the choice of parameters. Quadratic local interpolants might yield better results, but the method is already very 8 a) b) c) d) Figure 3: Deformation with inverse distance weighted methods. a) Original Shepard approach without local interpolant. b) Linear local interpolant, = 2. c) Linear local interpolant, = 4. d) Locally bounded weight functions, R = 256. time{consuming with linear interpolants. With quadratic interpolants, this drawback would become more severe. Thus, even though the method is exible, it does not appear very suitable for image warping. 9 5 Radial Basis Functions Another popular approach to scattered data interpolation is to construct the interpolant as a linear combination of basis functions and then to determine the coe cients of the basis functions, f(x) = n X i=1 i R(di(x)) + pm(x): (11) The values of the basis function R depend only on the distance from the data point, and are thus called radial. pm(x) is a polynomial of degree m. It assures a certain degree m of polynomial precision. The coe cients i are calculated by putting the data points into (11) and solving the resulting system of linear equations. The di erentiability of this interpolation method depends directly on the di erentiability of the basis functions R used. Well{known radial basis functions are multiquadrics, originally proposed by R. Hardy [10, 11], R(d) = (d2 + r2) =2 with r > 0 and 6= 0: (12) Hardy's original approach does not contain a polynomial pm. Therefore, it has a polynomial precision of 0. For the exponent , Hardy proposes = 1. = 1 has also been used successfully. For r 6= 0, the basis functions are in nitely di erentiable, so the resulting interpolation is also in C1. The characteristic radius r > 0 can be chosen arbitrarily. It determines the smoothness of interpolation at the given data points. For image deformation, the selection of r is critical for good results. Values too small lead to undesirable unevenness in the deformed image, values too big can lead to foldover. With xed values of r, both e ects are likely to occur even within the same image, as shown in g. 4a, where foldover is evident near the eye and, at the same time, the rim of the hat is unpleasantly bent. Thus, it appears necessary to choose individual values of r for each control point. Following a suggestion in [2], we use individual values ri for each data point pi, computed from the distance to the nearest neighbor, i.e. ri = mini6=j di(xj): (13) 10 This causes the deformation to be softer when data points are widely spaced and stronger when they are closer together. The application of radial basis functions to the problem of deformation gives f (p) = n X i=1 i R(di(p)) + pm(p): (14) Now i 2 IR2 are coe cient vectors and pm : IR2 ! IR2 is a function with polynomial components of degree m. It turns out that Hardy's multiquadrics approach without a polynomial term pm does not yield satisfactory results as it produces distortions even if all control points remain at their original positions and diverges quickly outside the convex hull of the control points. This behavior is caused by the non{linearity of the basis functions. Introducing a linear pm, i.e. m = 1, as a global transformation with coe cients computed by the least squares method to minimize errors for the known control points solves this problem. More easily, under most circumstances, an identical transform is su cient. The resulting mapping function is then f (p) = n X i=1 i ((di(p))2 + r2 i ) =2 + p (15) The image in g. 4d has been computed with this method, with = 1. Fig. 4c shows a test grid deformed with the same parameters to demonstrate the smoothness of the deformation. A disadvantage of radial basis functions { like all global interpolation methods { is that for each pixel, all control points have to be taken into account. Thus, the algorithmic complexity is O(n N ), where n is the number of control points and N the number of pixels to be displaced, as opposed to triangulation based methods with a complexity O(N ). There is, however, room for improvements in computational e ciency in radial basis functions. For radial basis functions with = 1, n N square roots have to be computed. This is undesirable even if a fast square root approximation is used. Fortunately, our experiments show that = 2 gives results which are almost as good, as shown in g. 4b. With this exponent, the square root is replaced by a division. Additionally, if the pixel locations of the input image are on a regular grid, and the initial locations of the control points also lie on that 11 grid, the number of distances that has to be evaluated can be greatly reduced by computing the pixel displacements per control point radially outward from the control point, so that the typically eight grid points with the same distance from the control point are treated simultaneously. The performance gain from this algorithm as compared to a straightforward scanline algorithm is, however, largely o set by the larger number of boundary checks necessary. a) b) c) d) Figure 4: Deformation with radial basis functions. a) Fixed r = 50, = 1. b) Adaptive r, = 2. c) Test grid deformed with adaptive r, = 1. d) Image deformed with the same parameters. 12 6 Locally Bounded Radial Basis Functions The algorithm outlined above, where evaluation of pixel displacements spreads outwards from the control points, also opens the door for another more signi cant improvement. As pointed out above, multiquadrics with a value of = 2 give quite reasonable results for image deformation. These vanish for large distances. It therefore seems plausible to limit the range of in uence of the basis functions, so that evaluation stops at some distance from the control points. This way it is possible to beat the O(n N ) complexity limit. The range of in uence can be limited by multiplication with a mollifying function, so that the radial basis function has the form Ri(d) = (d2+r2 i ) (1 (d=Li)2) + with x + = fx ; 0; if x > 0 if x 0 (16) where Li is the limit of the range of in uence of the control point pi. The resulting basis functions are in C 1. In our experiments, a simpler and more e cient approach has been used. We subtract a constant o set from the basis functions, i.e. Ri(d) = (d2 + r2 i ) i + (17) where i is chosen so that Ri(Li) = 0. The derivatives of this function are not de ned for d = Li, so that Ri 2 C0, but for su ciently large Li, this does not matter in practice. Two strategies for determining the Li have been tried. The rst strategy just sets Li to a multiple of ri, the characteristic radius of the multiquadrics. This ensures some overlap of the in uence ranges of near neighbors. However, it can lead to gaps where no deformation takes place if control points are spaced very irregularly. The second strategy sets Li to a distance where the in uence of Ri is below a certain threshold T . To achieve this, the coe cients i are rst evaluated for normal multiquadrics (eq. 12). The o set i is then set to i = T=max(j i;xj; j i;yj). This way, the in uence of the modi ed radial basis function vanishes where the original multiquadric would have had an in uence of at most T . The coe cients i then have to be re{evaluated for the modi ed multiquadrics with the o set i. This could be iterated 13 a) b) Figure 5: Deformation with locally bounded radial basis functions. a) Locally bounded with Li 2ri and T = 5. b) Locally bounded with Li 1:4ri and T = 10. several times, but from our experience, one iteration step appears to be quite su cient. The second strategy has one major aw in that it is possible to get Li = 0, which would cause the linear equation system for the i to have no solution. This problem is solved by specifying a minimum range Lmin so that Li Lmin. For Lmin, we use a value determined with the rst strategy. This combination of the two strategies also helps to obtain greater stability with regards to unwanted artefacts in the deformation, like foldover or local roughness. Our experiments show that with locally bounded radial basis functions, speedups in the order of 10 can be achieved without serious loss of quality. Fig. 5 shows examples of locally bounded multiquadrics. Fig. 5a shows the result with a threshold of T = 5 pixels and a minimum range of 2ri. There is no visible di erence to multiquadrics with = 2, but the computation was about seven times faster. Fig. 5b shows deformation with T = 10 and a minimum range of 1:4ri. Despite the fact that the time required for the deformation is halved again, the quality is quite striking. The only apparent problem is the slight foldover near the tip of the hat. 14 7 Performance Comparisons Table 1 shows the result of some run{time measurements with the example picture used throughout this article. This is a 512 512 pixel, 256 gray levels image, deformed with 29 control points. The measurements took place on a Sun Sparc2. The measured time value is the user time spent on behalf of the deformation, as determined from calls to the system function rusage before and after the call to the deformation routine. This time measurement is relatively una ected by system load. No le operations take place during the measurement. Detailed pro ling shows that very little of the time needed for deformation is spent for preparational operations like solving the equation system for multiquadrics. Almost all the time is used for the actual evaluation of the transformed locations of the pixels. deformation method average time triangulation, corners 0.9 sec triangulation, poles 1.6 sec inverse distance, = 2 75.6 sec multiquadrics, scanline, = 1 92.4 sec multiquadrics, scanline, = 2 29.4 sec multiquadrics, radial, = 2 24.8 sec multiquadrics, bound = 3 4.9 sec multiquadrics, bound = 2, thresh = 5 3.4 sec multiquadrics, bound = 1.4, thresh = 10 1.7 sec resampling with bilinear ll 32.1 sec Table 1: Run{time measurements for the various algorithms. For triangulation based methods, four more control points have to be added to ensure that the convex hull of the control points covers the whole image. In the \corners" measurement, these points were in the corners of the image. In the \poles" measurement, the points were positioned far away in the north, east, west, and south directions of the image. This ensures that the interior triangulation is not a ected by the added points, but adds some overhead for points outside the image. This second set of additional points was used for the images shown. 15 For inverse distance weighted methods with = 2, no exponentiationis necessary. Other exponents give lower performance.For multiquadrics with = 1, the square root function from the stan-dard math library is used. For, = 2, no exponentiation is necessary.The result of the deformation with locally bounded radial basis functionswith xed Li = 3 ri is very similar to that shown in g. 5a, but takes moretime to compute.There is a number of possibilities to further increase the speed of com-putation. One approach would be to reduce the number of points mappedwith a global method by skipping pixels and mapping the pixels in be-tween with a simpler method, e.g. a bilinear interpolation.Also, if displacements are computed with a scanline algorithm, theoperation to be executed is the same for all pixels, so that global interpo-lation methods, and radial basis functions in particular, are ideally suitedfor parallelization.In an interactive environment, locally bounded multiquadrics couldbe further modi ed so that when a new control point is added, only thecoe cient for this new control point is evaluated to displace this pointas desired, and the complete system of equations re{evaluated only whentime permits. The displacements caused by a single point can be evaluatedin a small fraction of a second.However, with the current implementation of the resampling process,this by far dominates the overall time needed for warping, so that furtherimprovements of the performance of the deformation stage do not appearuseful at the moment. The time needed for resampling the test image isgiven in the last row of table 1.8 ConclusionMapping functions for image warping based on scattered data interpola-tion methods have been described. Methods based on radial basis func-tions, particularly multiquadrics, have been shown to be suitable for ap-plication to image warping. These functions produce deformations in C1,and are relatively insensitive to foldover problems. For strong deforma-tions, where foldover is more likely to occur, the method lends itself wellto decomposition into step{wise deformation [17].16 A disadvantage of global deformation methods relative to the morecommon triangulation based methods is the higher algorithmical com-plexity. A local method based on multiquadrics has been presented thatovercomes this problem without sacri cing quality. The performance ofthis method can compete with that of triangulation based methods and islikely to surpass them if smooth interpolation is used within the triangles.Further work will be devoted to the application of these methods totwo{dimensional morphing and to the interpolation between planar slicesin medical data sets. Also, since the methods described here depend onlyon Euclidean distance, they can easily be adapted to volume deformationapplications, e.g. the registration of voxel data.References[1] P. Buche and J. Camillerapp. Serial cuttings matching: An appli-cation to muscle ber characterization. In Proc. EUROGRAPHICS'91, pages 329{340, 1991.[2] M. Eck. Interpolationsmethoden zur Rekonstruktion von 3D{O-ber achen aus ebenen Schnittfolgen. CAD und Computergraphik,13(5):109{120, Feb. 1991.[3] R. Farwig. Rate of convergence of Shepard's global interpolationformula. Mathematics of Computation, 46(174):577{590, 1986.[4] R. Franke and G. Nielsen. Smooth interpolation of large sets ofscattered data. Int. Journal for Numerical Methods in Engineering,15:1691{1704, 1980.[5] R.C. Gonzalez and P. Wintz. Digital Image Processing. AddisonWesley, Reading, MA, 2nd edition, 1987.[6] W.J. Gordon and J.A. Wixom. Shepard's method of \metric inter-polation" to bivariate and multivariate interpolation. Mathematicsof Computation, 32(141):253{264, 1978.[7] A. Goshtasby. Piecewise linear mapping functions for image registra-tion. Pattern Recognition, 19(6):459{466, 1986.17 [8] A. Goshtasby. Piecewise cubic mapping functions for image registra-tion. Pattern Recognition, 20(5):525{533, 1987.[9] A. Goshtasby, G.C. Stockman, and C.V. Page. A region{based ap-proach to control point selection with subpixel accuracy. IEEE Trans.Geosci. Remote Sens., GE-24:390{399, 1986.[10] R.L. Hardy. Multiquadric equations of topography and other irregu-lar surfaces. J. Geophys. Res., 76:1905{1915, 1971.[11] R.L. Hardy. Theory and applications of the multiquadric{biharmonicmethod. Computers and Mathematics with Applications, 19:163{208,1990.[12] P.S. Heckbert. Survey of texture mapping. IEEE Computer Graphicsand Applications, 6(11):56{67, 1986.[13] L.N. Kanal, B.A. Lambird, D. Levine, and G.C. Stockman. Digitalregistration of images from similar and dissimilar sensors. In Proc.Int. Conf. Cybern. Society, pages 347{351, 1981.[14] P. Lancaster and K. Salkauskas. Surfaces generated by moving leastsquares methods. Mathematics of Computation, 37(155):141{158,1981.[15] C.L. Lawson. Software for C1 surface interpolation. In J.R. Rice,editor, Mathematical Software III, pages 161{194. Academic Press,1977.[16] J. Rasure and M. Young. An open environment for image processingsoftware development. In 1992 SPIE/IS&T; Symposium on ElectronicImaging, SPIE Proceedings Vol. 1659, February 1992.[17] D. Ruprecht and H. Muller. Free form deformation with scattereddata interpolation methods. Universitat Freiburg, Institut fur Infor-matik, Bericht 41, 1991. To be published in Computing Suppl. 8,Springer{Verlag, Wien 1992.[18] D. Shepard. A two{dimensional interpolation function for irregularlyspaced data. In Proc. 23 Nat. Conf. ACM, pages 517{524, 1968.18 [19] A.R. Smith. Planar 2-pass texture mapping and warping. ComputerGraphics, 21(4):263{272, 1987.[20] G.C. Stockman, S. Kopstein, and S. Benett. Matching images tomodels for registration and object detection via clustering. IEEETrans. Pattern Anal. Mach. Intell., 4:229{241, 1982.[21] G. Wolberg. Digital Image Warping. IEEE Computer Society Press,1990.[22] G. Wolberg and T.E. Boult. Separable image warping with spatiallookup tables. Computer Graphics, 23(3):369{378, 1989.